2-2

Physicist's Introduction to Mathematica 6

Third Edition 2007

PHYSICS 200

Nicholas Wheeler

Laboratory 3

Some Useful Standard Packages

The power of Mathematica is greatly enhanced by the large number of Standard Packages that are tucked into its memory, and stand available for activation on an as-needed basis. To find out about those, go to the Documentation Center and search for Standard Packages, which provides links to sources of more specialized information.  The Mathematica Packages provides general information. Standard Package Compatibility Guide opens to reveal a set of eleven topical heads, each of which opens to provide links to a large collection of sites, each of which provides updated information about a specific package.  Many other "non-standard packages" (some proprietary, some not) have been developed to serve specialized needs, and can be discovered on the web: ask Google. Here I can attempt only to illustrate the kind of material that lies hidden away in packages of various descriptions.

Physical Constants

Go to the Documentation Center and search for PhysicalConstants/tutorial/PhysicalConstants. To install the package, command

Needs["PhysicalConstants`"]

You now have at your command such information as the following:

PlanckConstant

SpeedOfLight

EarthRadius

FineStructureConstant

The fine structure constant ((e^2)/ℏc) is dimensionless (has the same value whatever units are employed), and has approximately the value 1/137:

N[(1/137)]

Units Conversions

To enable conversion  from one system of units to another on has available the resources of the Units Package. To learn about it, search for Units/tutorial/Units. To activate it, command

Needs["Units`"]

This puts you in position to enter such commands as the following:

Convert[95 Calorie, Erg]

Convert[16 Mile,Meter]

N[%]

Convert[1Year,Second]

How accurate is the claim that "one year is about π (10^7) long"?

(31536000/(π (10^7)))//N

Convert[ProtonMass,AtomicMassUnit]
Convert[NeutronMass,AtomicMassUnit]

Roughly how many such particles would it take to fabricate the Sun?

Convert[SolarMass,AtomicMassUnit]

MKS[65Mile/Hour]

N[%]

CGS[65 Mile/Hour]//N

Calendar

To read about this surprisingly useful "Miscellaneous Standard Package," launch a search for Calendar.  To activate the package, command

Needs["Calendar`"]

and I now provide quick illustration of its use. My friend Oya was born on 15 March 1944.

DayOfWeek[{1944,3,15}]

we learn that she was born on a Wednesday. I am writing on 31 July 2007, and from

DaysBetween[{1944,3,15},{2007,7,31}]

we learn that she is today 23,148 days old. That comes to

Convert[23148 Day,Second]

so she will—early tomorrow—be able to celebrate a moment when she has become precisely 2 billions seconds old:

Convert[2000000000Second,Day]//N

Elementary Algebra

Short List of Basic Operations

Mathematica has lots to say about algebraic subjects of all descriptions, many quite esoteric. My simple intent here is to discuss the resources that come into view when one goes to Palettes > AlgebraicManipulation, which can be read about by consulting on or another of the eleven tutorial notebooks to which tutorial/AlgebraicManipulationOverview provides links.

Mathematica looks upon (((x+y))^7)as a primitive object.

(((x+y))^7)

But it will expand that binomial if you ask it to:

Expand[(((x+y))^7)]
Expand[(2+3x+17(x^2))((x^2)-9)((x^2)+25)]
Expand[(2+3x+17(x^2))((x^2)-9)((x^2)+25)]//TraditionalForm

Standardly, the powers (x^n) are presented in ascending order.  The command //TraditionalForm reverses the order and adjusts the font.

The action of Together is similarly self-explanatory, and opposite to that of Apart:

Clear[a,b,c,d]

Together[(a/b)+(c/d)]

Apart[((b c+a d)/(b d))]

REMARK: It is often more convenient to use post-commands to achieve the same effect:

(a/b)+(c/d)//Together

%//Apart

The Apart command is used to produce the partial fraction representation of ratios, which is often useful in pure/applied mathematical work:

Apart[(1/((x-a)(x-b)))]

One uses partial fractions to "separate the singularities," as in the following little example:

Plot[(1/((x-1)(x-2))),{x,0,3},PlotRange->{-40,40}, Ticks->{{1,2,3},{-40,-20,20,40}}]

(1/((x-1)(x-2)))//Apart

Plot[{(1/(-2+x)),-(1/(-1+x))},{x,0,3},PlotRange->{-40,40}, Ticks->{{1,2,3},{-40,-20,20,40}},PlotStyle->{Red,Blue}]

Three very useful commands are

?Cancel
?ExpandDenominator
?ExpandNumerator

—the use of which I now illustrate:

(((1+x)(2+x)(7+x))/((1+x)(3+x)(1+2x+3(x^2))))//Cancel

(((1+x)(2+x)(7+x))/((1+x)(3+x)(1+2x+3(x^2))))//ExpandDenominator

(((1+x)(2+x)(7+x))/((1+x)(3+x)(1+2x+3(x^2))))//ExpandNumerator

One should not confuse the number-theoretic command

?FactorInteger

which works this way

FactorInteger[284]

(2^2)*71

…with the command

?Factor

which possesses several optional variants

Options[Factor]

and works this way:

Factor[(x^2)-9]

Factor[(x^2)+9]
Factor[(x^2)+9, GaussianIntegers->True]

But factoring 2+3x+17(x^2) leads away from the (Gaussian) integers—amounts to looking for roots of the polynomial:

Finding the Roots of a Polynomial

The celebrated Fundamental Theorem of Algebra (Gauss) asserts that every polynomial of degree n (coefficients real or complex) possesses n roots (real or complex).  Those can be discovered by means of the following commands:

?Solve
?NSolve

Solve[2+3 x+17 (x^2)==0,x]

NSolve[2+3 x+17 (x^2)==0,x]

The quadratic formula has been known from dim antiquity

Solve[a (x^2)+b x+c==0,x]

Given a polynomial of order n, one can without loss of generality always arrange to have the the coefficient of (x^n) be unity, and the coefficient of (x^(n-1)) be zero. Assume those adjustments to have been made.  In the 16th Century Jerome Cardan, Niccolo Tartaglia and Ludovico Ferrari carried the problem of solving the cubic and quartic as far as one could without knowledge of complex numbers (the invention of which was stimulated by their success):

Solve[ (x^3)+b x+c==0,x]

Solve[ (x^4)+a (x^2)+b x+c==0,x]

All efforts to solve the quintic (which everyone expected to be difficult) failed…for what Galois and Abel showed (early 19th Century) to be this very good reason: solution (in the sense of the preceding formulae) is impossible, except in favorable special cases. Mathematica has this to say in the general case:

Solve[ (x^5)+w (x^3)+a (x^2)+b x+c==0,x]

But (approximate) numerical solution is possible in all specific cases:  the preceding command supplies an uninformative result

Solve[(x^5)+ (x^3)+2 (x^2)+3 x+4==0,x]//TableForm

but look to its numerical sybling:

NSolve[(x^5)+ (x^3)+2 (x^2)+3 x+4==0,x]//TableForm

({{x->-1.1181754240534836`}, {x->-0.4326913321786198`-1.0949614143362267` I}, {x->-0.4326913321786198`+1.0949614143362267` I}, {x->0.9917790442053616`-1.2637504112999582` I}, {x->0.9917790442053616`+1.2637504112999582` I}})

Here is a more legible version of the same statement:

NSolve[(x^5)+ (x^3)+2 (x^2)+3 x+4==0,x,6]//TableForm

I introduce now a graphic function, the design of which was adapted from Chapter 36 of Glynn & Gray's Beginner's Guide to Mathematica 4:

RootPlot[poly_,s_]:=ListPlot[Map[{Re[#1],Im[#1]} &,x/.N[Solve[poly==0,x]]], PlotRange->{{-s,s},{-s,s}}, AspectRatio->Automatic,PlotStyle->PointSize[.02],AxesLabel->{"Re","Im"}]

REMARK: It is surprising to me that in 1992, in an article entitled "Some elements of Mathematica design" (visit MouseOver) Wolfram himself wrote at length about why he had avoided the introduction of such an obviously useful command.

Here poly refers to the polynomial whose roots one wants to display as points on the complex plane, and s sets the scale of the figure (PlotRange). The command works like this:

RootPlot[(x^5)+ (x^3)+2 (x^2)+3 x+4,2]

Another Example

Let

p[x_]:=(x^6)-2(x^2)+4x+4

and look to the roots of p[x]=0:

RootPlot[p[x],1.5]

When we plot the polynomial

Plot[p[x],{x,-1.5,1.5}]

we see the real roots as axis-crossings, but are provided with only an obscure hint concerning the location of the complex roots.

Example: The 17th Roots of Unity

RootPlot[(x^17)-1,1.5]

Quite a similar figure is seen to result from

RootPlot[1+x+(x^2)+(x^3)+(x^4)+(x^5)+(x^6)+(x^7)+(x^8)+(x^9)+(x^10)+(x^11)+(x^12)+(x^13)+(x^14)+(x^15)+(x^16),1.5]

The graphs are so similar because the two polynomials are very closely related:

Series[(((x^17)-1)/(x-1)),{x,0,17}]

Removal of the (x-1) factor has erased one of the dots (roots). The remaining 16 roots are the 17th roots of unity except for unity itself:

Solving Systems of Simultaneous Equations

Generic Pair of Inhomogeneous Linear Equations

Clear[a,b,c,d,p,q,x,y]

Solve[{a x+b y==p, c x+d y==q},{x,y}]

Here we ask for an alternative display of the same output:

Solve[{a x+b y==p, c x+d y==q},{x,y}]//TableForm

Triple of Inhomogeneous Linear Equations with Numeric Coefficients

Solve[{1 x+ 2 y+3z==p,4 x+5 y+6z==q,7x+8y+9z==r}, {x,y,r}]//TableForm

A Nonlinear Coupled Pair

Solve[{(x^2)+(y^2)==2,x+y==1},{x,y}]//TableForm

To obtain the numerical evaluation of that result, we might enter the post-command

N[%]

or might back up and modify the initial command:

NSolve[{(x^2)+(y^2)==2,x+y==1},{x,y}]//TableForm

Here is a graphic representation of what Mathematica has just accomplished:

ContourPlot[{(x^2)+(y^2)==2,x+y==1},{x,-2,2},{y,-2,2}, Frame->False, Axes->True]

A Case in which "No Solution Exists"

ContourPlot[{(x^2)+(y^2)==1/3,x+y==1},{x,-2,2},{y,-2,2}, Frame->False, Axes->True]

It is evident to the eye that in this case the curves do not intersect = the system "has no solution." But the command

NSolve[{(x^2)+(y^2)==1/3,x+y==1},{x,y}]//TableForm

reminds us that what we really mean to say is that the system "has no real solution."

Application to the Theory of One-dimensional Elastic Collisions

The pre-collisional velocities of M and m are U and u. If the collision is elastic, then both energy and momentum are conserved:  to discover the post-collisional velocities V and v we must solve the coupled system

M (V^2)+m (v^2)=M (U^2)+m (u^2)

M V+m v=M U+m u

This (on account of the nonlinearity) is a bit tricky when approached analytically, but is duck soup for Mathematica:

Solve[{M (V^2)+m (v^2)==M (U^2)+m (u^2), M V+m v==M U+m u}, {V, v}]//TableForm

The first solution is trivial/spurious, the second describes the effect of a collision, and leads us to define "post-collisional velocity functions"

V[U_,u_,M_,m_]:=((M-m)/(M+m))U+((2m)/(M+m))u
v[U_,u_,M_,m_]:=((2M)/(M+m))U-((M-m)/(M+m))u

First illustrtative application: Suppose that M = m and that U = 0.  We then have

List[V[0,u,m,m],v[0,u,m,m]]

according to which the identical particles have simply traded velocities.

Second illustrative application: Suppose now that M = 1000 m and that (as before) U = 0. We get

List[V[0,u,1000m,m],v[0,u,1000m,m]]

The massive particle (initially at rest) is set into slow motion, while the incident particle rebounds with essentially the negative of its initial velocity.

We look now to an adjustable graphical representation of the algebraic situation.

Manipulate[ContourPlot[{3(x^2)+(y^2)==2w,3x+y==p},{x,-2,2},{y,-2,2}, Frame->False, Axes->True],{{w,1,"energy"},.2,2},{{p,1,"momentum"},-3,3}]

It is interesting to note that in cases where line is tangent to ellipse, the point of tangency always occurs on the 45º "equi-velocity" line, for this very good reason: if particles are moving with the same velocity then a "collision" is an inconsequential event, a lingering kiss.

If the collision were inelastic then, though the momentum line would remain invariant, the energy ellipse would shrink to tangency, by which time the system would have given up (thermalized) the maximal allowed amount of kinetic energy.

Trigonometry

Looking again to the AlgebraicManipulation palette, we are reminded that commands like Expand[■], Factor[■] and Simplify[■]are of a type known as formula manipulation commands. We experiment now with some trigonometrically informed relatives of those. The simple commands in question are

?TrigExpand
?TrigFactor
?TrigReduce
?TrigToExp
?ExpToTrig

A few examples will make clear what they do:

TrigExpand[Cos[α+β]]

TrigFactor[Cos[α]+Cos[β]]

TrigExpand[%]

Simplify[%]

TrigReduce[((Cos[α])^3)]
TrigReduce[Cos[α]Sin[β]]

TrigFactor[%]

TrigToExp[((Sin[α])^3)]

ExpToTrig[(E^(I θ))]

Simplify[((Sin[α])^2)+((Cos[α])^2)]

With the aid of the preceding commands one can—with experience—reconstruct and enlarge upon the standard lists of trigonometric identities.  Experience is required because any given trigonometric expression can be written in a zillion ways, and one must learn how to tell Mathematica the goal one has in mind.

And in the latter regard one must be braced for some surprises.  For example: suppose that—with

TrigExpand[(((Sin[α]+Sin[β]))^2)]

in mind—we were to ask for

TrigExpand[(((Sin[10 t]+Sin[11 t]))^2)]

This is an unanticipated result.  What's going on?  Mathematica interprets  Sin[10 t] to mean Sin[t+t+t+t+t+t+t+t+t+t] because 10 is an integer.  To get the result we anticipated we must ask Mathematica to think of 10 as just a real number:

TrigExpand[(((Sin[10.0 t]+Sin[11.0 t]))^2)]

Matrices

For introductory information, go to Documentation Center, launch a search for guide/MatrixOperations and examine the TUTORIALS listed at the bottom of the page.  A matrix is a list of lists, so

{{a,b},{c,d}}
{a,b,c,d}
{{a},{b},{c},{d}}

are all matrices. To produce objects that look more like matrices, enter the commands

MatrixForm[{{a,b},{c,d}}]
{{a,b},{c,d}}//MatrixForm

{{a,b,c,d}}//MatrixForm
{{a},{b},{c},{d}}//MatrixForm

Creating Matrices

The simplest way to make a matrix appear at the cursor is to proceed Insert > Table/Matrix and select New. [REMARK: In Mathematica 5 one could achieve the same result from the keyboard, by typing ⇧⌘C, but in Mathematica 6 that capability has been lost.] A menu will appear, where you will be asked to state what it is you desire (select amongst several options). Do so, then hit OK. You will have created something like this:

(({{\[Placeholder], \[Placeholder], \[Placeholder], \[Placeholder]}, {\[Placeholder], \[Placeholder], \[Placeholder], \[Placeholder]}, {\[Placeholder], \[Placeholder], \[Placeholder], \[Placeholder]}, {\[Placeholder], \[Placeholder], \[Placeholder], \[Placeholder]}}))

Place the cursor on the 11 box, type a, then type ⇥ b ⇥ c ⇥ d… to produce the following:

Clear[a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p]

(({{a, b, c, d}, {e, f, g, h}, {i, j, k, l}, {m, n, o, p}}))

The column-adding effect of placing the cursor between a and b and typing Ctrl+comma is shown  below

(({{a, \[Placeholder], b, c, d}, {e, \[Placeholder], f, g, h}, {i, \[Placeholder], j, k, l}, {m, \[Placeholder], n, o, p}}))

while the following matrix shows the row-adding effect of placing the cursor as before, but typing Ctrl+⏎

(({{a, b, c, d}, {\[Placeholder], \[Placeholder], \[Placeholder], \[Placeholder]}, {e, f, g, h}, {i, j, k, l}, {m, n, o, p}}))

To remove/modify an element, select that element and hit ⌫; thus

(({{a, b, c, d}, {e, f, \[Placeholder], h}, {i, j, k, l}, {m, n, o, p}}))

To remove a row (or column), select that row/column and hit ⌫; thus

(({{a, c, d}, {e, g, h}, {i, k, l}, {m, o, p}}))

TIP: It rapidly becomes tedious to return again and again to the Table/Matrix menu to create successive matrices of similar structure. In preparing the preceding paragraph, I proceeded in the manner just described to create

(({{□, □, □, □}, {□, □, □, □}, {□, □, □, □}, {□, □, □, □}}))—actually, to create (({{a, b, c, d}, {e, f, g, h}, {i, j, k, l}, {m, n, o, p}}))—and then used Copy/Paste.

Notice that options available on the Table/Matrix menu permit one to create matrices with (for example) the structure

(({{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}))

with a single keystroke. Alternatively, one might command

IdentityMatrix[4]//MatrixForm

REMARK: When using ⇥ to fill in the elements of  (({{□, □, □, □}, {□, □, □, □}, {□, □, □, □}, {□, □, □, □}})) it is useful to remember that the control keys ↑, ↓, → and ← act to relocate the "selected box" in obvious ways.

Matrices can constructed at the keyboard, and this is sometimes the fastest way to go. An example serves to describe the simple routing. In Input mode, type ()← to create

()  with the cursor blinking inside

Now type aCtrl,bCtrl,cCtrl,d (NOTE: Ctrl, means "type comma with the Ctrl key depressed") to produce

(({{a, b, c, d}}))  with the cursor blinking after the d

Now type Ctrl⏎ to produce

(({{a, b, c, d}, {\[Placeholder], \[Placeholder], \[Placeholder], \[Placeholder]}}))

Fill in the new row. How rows can be lengthened/shortened, new rows can be created and old ones elminated should by now be clear.

Subscripts & Superscripts/Exponents

Subscripts are useful in many contexts, but in matrix algebra they are especially important, since they serve very efficiently  to make it clear who stands where, and to clarify the relatinships that develop. For information, search for tutorial/EnteringTwoDimensionalExpressionsOverview and select Typing Superscripts.

To create

(Subscript[a, 1])

one might type a and then hit the (Subscript[□, □]) button near the bottom of the BasicMathInput palette.  I find it quicker, however, to type aCtrl _1Ctrl⇧, where the final Ctrl⇧ is needed  to get the cursor out of the basement.

Double superscripts are a bit more ticklish. To create

(Subscript[a, 12])

you might proceed essentially as before: type aCtrl _12Ctrl⇧. And for many purposes that is sufficient, as I demonstrate:

A=(({{(Subscript[a, 11]), (Subscript[a, 12])}, {(Subscript[a, 21]), (Subscript[a, 22])}}));
B=(({{(Subscript[b, 11]), (Subscript[b, 12])}, {(Subscript[b, 21]), (Subscript[b, 22])}}));

A
B

NOTE that Mathematica insists upon thinking of matrices as lists of lists. We look to the matrix product of A and B:

A.B//MatrixForm

Contexts arise, however, in which one wants Mathematica to think of the subscript on

(Subscript[a, 12])

not as a unitary symbol ("twelve" in this instance) but as "one, two". One has then to separate the "one" from the "two" with a comma:

(Subscript[a, 1,2])

If one finds that too distracting in a context like

(({{(Subscript[a, 1,1]), (Subscript[a, 1,2])}, {(Subscript[a, 2,1]), (Subscript[a, 2,2])}}));

%//MatrixForm

one can make use of the invisible comma, which is created Esc,Esc. One then has

(({{(Subscript[a, 1,1]), (Subscript[a, 1,2])}, {(Subscript[a, 2,1]), (Subscript[a, 2,2])}}))//MatrixForm

Which looked fine as input, but which as output reacquired its former ugliness.

One might suppose that spaces would do the job, but no:

(({{(Subscript[a, 1 1]), (Subscript[a, 1 2])}, {(Subscript[a, 2 1]), (Subscript[a, 2 2])}}))//MatrixForm

Mathematica interprets xy to be a unitary symbol, but x y as a command to construct the product of x and y:

x y===x*y

IMPORTANT NOTE: In this respect, matrices are exceptional. To signal that matrices A and B are to be multiplied, neither A B  nor A*B works. One must write A.B (as above).

Superscripts

Identical remarks pertain the creation of superscripted matrices, the only difference being that in place of  aCtrl _12Ctrl⇧ we type pCtrl ^12Ctrl⇧. We have

P=(({{(p^11), (p^12)}, {(p^21), (p^22)}}));
Q=(({{(q^11), (q^12)}, {(q^21), (q^22)}}));

P
Q

P.Q//MatrixForm

No surprises. But consider:

A.A//MatrixForm
P.P//MatrixForm

Mathematica has interpreted the superscripts to be exponents, and has combined them by the Law of Exponents!!

MORAL: Superscripts are best avoided, not only in matrix algebra but in all algebraic contexts, for they occupy a space reserved for exponents.

Some Common Matrix Operations

Let the matrix A be defined as before:

A//MatrixForm

Then to evaluate the trace we command

Tr[A]

To evaluate the determinant we command

Det[A]

To construct the transpose we command

Transpose[A]//MatrixForm

To construct the inverse we command

Inverse[A]

and check that this matrix has the required property:

Inverse[A].A//Simplify//MatrixForm

Let us now define

rowvector=(({{(Subscript[a, 1]), (Subscript[a, 2]), (Subscript[a, 3])}}))
columnvector=(({{(Subscript[b, 1])}, {(Subscript[b, 2])}, {(Subscript[b, 3])}}))
squarematrix=(({{(Subscript[m, 11]), (Subscript[m, 12]), (Subscript[m, 13])}, {(Subscript[m, 21]), (Subscript[m, 22]), (Subscript[m, 23])}, {(Subscript[m, 31]), (Subscript[m, 32]), (Subscript[m, 33])}}))

The  following command is only occasionally useful, for it responds to a question that is seldom in doubt:

Dimensions[rowvector]
Dimensions[columnvector]
Dimensions[squarematrix]

The matrix product

rowvector.columnvector

%//MatrixForm

is a 1×1 matrix, a list with a single element, that element being itself a list with a single element, that element being familiar as the dot product. We can use a standard list manipulation technique to peal alway the distracting brackets:

%[[1]]

%[[1]]

More compactly:

(rowvector.columnvector)[[1]][[1]]

REMARK: Just about everything, if you penetrate deeply enough into the bowels of Mathematica, becomes a "list" of one sort or another, and elaborate means are provided for the manipulation of lists. Sooner or later you will need to acquire familiarity with some of that material. For basic information, search for tutorial/ListsOverview and look at the links provided there.

Reversing the order of the matrix product studied above, you obtain

columnvector.rowvector//MatrixForm

which is a 3×3 matrix. Matrices of this specialized design are sometimes called "dyadics."

Here are a couple of matrix operations that—Mathematica is quick to inform you—make no sense:

rowvector.rowvector
columnvector.columnvector

Because matrix multiplication is in some respects a bit quirky, one must take special steps to raise a square matrix to a power. Comparing

(({{1, 2}, {3, 4}})).(({{1, 2}, {3, 4}}))//MatrixForm

with

(((({{1, 2}, {3, 4}})))^2)//MatrixForm

shows clearly that commands of the latter sort do not do the job (what we have in fact demonstrated is that squaring is a "Listable" operation: the meaning of "Listable" will be discussed in a moment). Available for such purposes is the command

?MatrixPower

Illustrative Application: Rotations in 2-Space

Define

S=(({{0, -1}, {1, 0}}));

Compare

S.S//MatrixForm
MatrixPower[S,2]//MatrixForm

Construct

R=(∑_(n=0)^6)(1/(n!))(θ^n)MatrixPower[S,n]//MatrixForm

Pretty clearly, the infinite sum gives

R=(({{Cos[θ], -Sin[θ]}, {Sin[θ], Cos[θ]}}));

which we recognize to be a rotation matrix:

Simplify[Transpose[R]==Inverse[R]]

Actually, Mathematica provides a command which eliminates the guesswork, and is very useful in a great variety of pure/applied contexts:

?MatrixExp

MatrixExp[θ S]//MatrixForm

Application of a Sort that Arises in Connection with the Solution of Differential Equations

Let

m[t_]:=Exp[a t]

Familiarly (if in a not-yet-familiar notation), m'[t] - a m[t] = 0:

m'[t] - a m[t] == 0

We look to the matrix analog of that elementary result: let

M[t_]:=MatrixExp[(({{-1.1, 0.9}, {-1.4, 0.3}}))t]

This is a fairly complicated matrix-valued function of t:

M[t]

But nevertheless,

M'[t]-(({{-1.1, 0.9}, {-1.4, 0.3}})).M[t]//Simplify//Chop//MatrixForm

where

?Chop

has served to remove computational artifacts of the order (10^(-17)).

Matrices as Arguments of "Listable" Commands

?Listable

Matrices are lists (of lists), so any command with the attribute "listable" will accept a matrix as its argument: thus

(({{1, x}, {(x^2), (x^3)}}))+a//MatrixForm
D[(({{1, x}, {(x^2), (x^3)}})),x]//MatrixForm
∫(({{1, x}, {(x^2), (x^3)}}))\[DifferentialD]x//MatrixForm

Characteristic Equation, Eigenvalues, Eigenvectors

I confine my initial illustrative remarks to the 2-dimensional case, and work with the generic matrix

M=(({{a, b}, {c, d}}));

The characteristic polynomial of the matrix is defined

Det[(({{a, b}, {c, d}}))-λ IdentityMatrix[2]]

but is achieved  much more simply the following command:

?CharacteristicPolynomial

CharacteristicPolynomial[M,λ]

The celebrated (and very often valuable) Cayley-Hamilton Theorem asserts—remarkably—that "every matrix satisfies its own characteristic equation." In the present instance we introduce the notations

unitmatrix=IdentityMatrix[2];
nullmatrix=(({{0, 0}, {0, 0}}));

and find that we do indeed have the matrix-valued identity

(-b c+a d)unitmatrix-a M-d M+M.M==nullmatrix

The "eigenvalues" of a square matrix are the roots of its characteristic polynomial. To discover those we might command

Solve[CharacteristicPolynomial[matrix,λ]==0,λ]

For the 2×2  M-matrix introduced just above this would give

Solve[CharacteristicPolynomial[M,λ]==0,λ]

but Mathematica provides a better way:

?Eigenvalues

Eigenvalues[M]//TableForm

Thus far we have asked Mathematica to do matrix-theoretic things we could as easily have done by hand. The power of Mathematica is more clearly evident in higher-dimensional situations. To have some representative material to work with, we construct a 5 by 5 matrix with random real-valued elements, making use of the command

?RandomReal

RandomReal[{-2,2},5]

NOTE: In a moment I am going to ask you to return to this point and to repeat the following sequence of commands.

★★★

randomrealmatrix=
{Table[RandomReal[{-2,2}],{5}],
Table[RandomReal[{-2,2}],{5}],
Table[RandomReal[{-2,2}],{5}],
Table[RandomReal[{-2,2}],{5}],
Table[RandomReal[{-2,2}],{5}]};

%//MatrixForm

To construct the determinant of such a matrix you would have to compute 5! five-number products and add the (appropriately signed) results; observe how quickly Mathematica does the job:

Det[randomrealmatrix]//Timing

Next we construct the characteristic polynomial

CharacteristicPolynomial[randomrealmatrix,λ]

The some instances, insignificantly small imaginary parts will appear in the coefficients, which are clearly computational artifacts. To get rid of them we would command

Chop[%]

Either by computation

Solve[%==0,λ]

or by means of the command

Eigenvalues[randomrealmatrix]//TableForm

we obtain the eigenvalues of our random 5×5 matrix. Notice that one (or sometimes three, sometimes five) is real, the others occur in conjugate pairs. Reactivating the definition (and making a notational adjustment λ → x as that definition requires)

RootPlot[poly_,s_]:=ListPlot[Map[{Re[#1],Im[#1]} &,x/.N[Solve[poly==0,x]]], PlotRange->{{-s,s},{-s,s}}, AspectRatio->Automatic,PlotStyle->PointSize[.02],AxesLabel->{"Re","Im"}]

we construct a graphic description of the spectrum of our 5×5 matrix; i.e., of its population of eigenvalues:

RootPlot[CharacteristicPolynomial[randomrealmatrix,x], 4]

Back up now to the point marked ★★★ and repeat that series of commands. Do that several times.  Observe that, while the numbers and final figure change randomly, the qualitative essentials of the situation persist.

Now construct the real symmetric matrix (the "symmetric part" of the current randomrealmatrix)

symmetricpart=(1/2)(randomrealmatrix+Transpose[randomrealmatrix]);

%//MatrixForm

The interesting thing about such matrices—revealed by the command

Eigenvalues[symmetricpart]//TableForm

is that their eigenvalues are all real:

RootPlot[CharacteristicPolynomial[symmetricpart,x], 4]

Real numbers that originate in this way—as the eigenvalues of a real symmetric matrix—are of fundamental importance to many branches of physics/engineering and pure/applied mathematics.

Now construct a 5 by 5 matrix with random complex-valued elements

✶✶✶

randomcomplexmatrix={
Table[RandomComplex[{-2-2I,2+2I}],{5}],
Table[RandomComplex[{-2-2I,2+2I}],{5}],Table[RandomComplex[{-2-2I,2+2I}],{5}],
Table[RandomComplex[{-2-2I,2+2I}],{5}],
Table[RandomComplex[{-2-2I,2+2I}],{5}]};

Construct also the hermitian part (complex generalization of the real symmetric part) of that matrix:

hermitianpart=(1/2)(randomcomplexmatrix+Transpose[Conjugate[randomcomplexmatrix]]);

%//MatrixForm

REMARK: You will have either to temporarily make the page wider or reduce the magnification of the display to see the whole matrix.

The calculation

Eigenvalues[randomcomplexmatrix]//TableForm

produces eigenvalues which now do NOT occur in conjugate pairs, while the numbers

Eigenvalues[hermitianpart]//Chop//TableForm

are, in every case, all real, as the following figure serves to dramatize:

RootPlot[CharacteristicPolynomial[hermitianpart,x], 5]

Return now to ✶✶✶ and repeat (several times) the preceding sequence of commands.

This result illustrates the general circumstance which, in fact and in particular, accounts for all the real numbers that issue theoretically from quantum mechanics, and that are supposed to account fo the results of "meter readings" (measurements).

The eigenvectors of a matrix arise by solution of an equation of the form

matrix.eigenvector==number*eigenvector

which is possible iff "number" is in fact one or another of the eigenvalues of the matrix. The eigenvectors can be produced by the command

?Eigenvectors

but it is usually more informative/efficient to use

?Eigensystem

The following simple example is designed to show how this works:

S=(({{1, 2}, {2, 3}}));

Eigensystem[S]

We use that information to define

eigenvalue1=2+(Sqrt[5])
eigenvector1=(({{-(3/2)+(1/2) (2+(Sqrt[5]))}, {1}}))
eigenvalue2=2-(Sqrt[5])
eigenvector2=(({{-(3/2)+(1/2) (2-(Sqrt[5]))}, {1}}))

and find that indeed

S.eigenvector1==eigenvalue1*eigenvector1
S.eigenvector2==eigenvalue2*eigenvector2

The evident orthogonality of the eigenvectors

Transpose[eigenvector1].eigenvector2//Simplify

Show[Graphics[{
{RGBColor[1,0,0],Arrow[{{0,0},Transpose[eigenvector1][[1]]}]},
{RGBColor[0,0,1],Arrow[{{0,0},Transpose[eigenvector2][[1]]}]}
}],AspectRatio->Automatic]

can readily be shown to be follow from the real symmetry of the matrix, just as the reality of the eigenvalues does.

Differential & Integral Calculus

For basic information concerning the calculus, in all of its many parts, as implemented in Mathematica, go to Documentation Center, search for guide/Calculus and NumericalCalculus/guide/NumericalCalculusPackage and inspect the links provided there. Here I will be able to discuss only a few of the most essential concepts.

Simple Indefinite/Definite Integrals

?Integrate

The ≫ takes you to ref/Integrate, which provides basic information, tutorials and links to sources that treat specialized topics.

Mathematica has its own esoteric ways of doing integrals, and has also committed to memory the contents of all standard integral tables. Integrals that it can't do analytically it is always prepared to do numerically.

Integral commands can be entered in three distinct ways: in "telegraphic style" from the keyboard; with 2D symbols constructed by keyboard codes; by means of the BasicMathInput palette. I will illustrate those three methods in turn.

Examples of Indefinite Integrals

Integrate[((Sech[x])^2),x]

That was a "telegraphic" command. We might alternatively have used the BasicMathInput palette (being very careful of its "unfortunate quirk") or typed EscinttEsc to produce

∫\[Placeholder]\[DifferentialD]\[Placeholder]

and then filled in the boxes to obtain

∫((Sech[x])^2)\[DifferentialD]x

Another indefinite integral :

∫(1/(1+(x^3)))\[DifferentialD]x

REMARK: To create fractions—like, for example,

(x/y)

—you can either use the BasicMathInput palette (being again very careful of its "unfortunate quirk") or—much more economically—type xCtrl/y, as instructed in tutorial/EnteringTwoDimensionalInput.

Examples of Definite Integrals

Look next to the related definite integral which I command first in telegraphic style

Integrate[1/(1+(x^3)),{x,0,1}]

and now in 2-dimensional style:

(∫_0^1)(1/(1+(x^3)))\[DifferentialD]x

REMARK: To create the box to receive the lower limit, type Ctrl- (control + hypen). For the upper limit type Ctrl%.

Another definite integral:

(∫_0^∞)6(E^(-a x))\[DifferentialD]x

NOTE: Don't fail to notice the 6 that stands n front of the If.

And yet another, the Gaussian integral:

(∫_(-∞)^∞)(E^(-a (x^2)))\[DifferentialD]x

It is one of Mathematica's more CURIOUS QUIRKs that in such output statements, the Assumptions stated are those which would cause the reported result to become invalid! Compare the results of the following commands, or which the first is quoted from Mathematica's output statement, and the second has the inequality adjusted to conform to our true intent and expectation:

Integrate[6(E^(-a x)),{x,0,∞},Assumptions->Re[a]≤0]
Integrate[6(E^(-a x)),{x,0,∞},Assumptions->Re[a]>0]

Here I use the Gaussian integral to provide the same sort of comparison:

Integrate[(E^(-a (x^2))),{x,-∞,∞},Assumptions->Re[a]≤0]
Integrate[(E^(-a (x^2))),{x,-∞,∞},Assumptions->Re[a]>0]

Functions Defined by Integrals

It is often the case even of integrals with simple integrands that they can be described only in terms of fancy functions (which in many case were historically defined by integrals): thus

∫((Log[1-(x^2)])/x)\[DifferentialD]x
∫Exp[1-(x^2)]\[DifferentialD]x
∫Sin[(x^2)]\[DifferentialD]x

But Mathematica stands ready to tell us about those higher functions: thus

?FresnelS

And often we can gain a sense of what they are about by constructing a figure:

Plot[FresnelS[x],{x,-10,10}]

Numerical Integration

But sometimes the arsenal of higher functions is inadequate to the task: some integrals are simply intractable. An example:

Integrate[(x^x),{x,0,1}]

In such cases one can, however, proceed numerically

NIntegrate[(x^x),{x,0,1}]

and can even plot the results of a continuum of numerical integrations:

Plot[NIntegrate[(x^x),{x,0,y}],{y,0,1}]

The numerical evaluation of integrals is, as we just saw, very easy to command, and Mathematica is able to deal automatically with various pathologies

Plot[(1/(Sqrt[x])), {x,0,1}]

NIntegrate[(1/(Sqrt[x])),{x,0,1}]

In some such cases Mathematica makes a best effort, but posts a cautionary flag:

Plot[(x^2)Sin[1/x], {x,-1,1}]
Plot[(x^2)Sin[1/x], {x,-0.1,0.1}]

NIntegrate[(x^2)Sin[1/x], {x,-1,1}]
NIntegrate[(x^2)Sin[1/x], {x,-1,2}]

For additional information about numerical integration, command

?NIntegrate

follow the » and inspect the links provided at the bottom of that page.

Ordinary & Partial Differentiation

?Derivative

The » takes you to ref/Derivative, which provides examples and useful links to tutorials and special topics.

Ordinary Derivatives (Functions of a Single Variable)

The preferred way to command low-order differentiation of a function of a single variable is to use primes:

f[x_]:=(x^n)

f'[x]
f''[x]
f'''[x]

We might alternatively have commanded

D[f[x],x]
D[f[x],{x,2}]
D[f[x],{x,3}]

although…

Partial Derivatives (Functions of Several Variables)

… D[ ] is more commonly used to construct ordinary derivatives of high order

D[f[x],{x,17}]

and partial derivatives:

g[x_,y_]:=(x^n) Sin[y]

D[g[x,y],x]
D[g[x,y],y]

D[g[x,y],{x,2}]
D[g[x,y],x,y]
D[g[x,y],{y,2}]

Expressions like D[h[x,y],x,x,x,x,y,y,y,y,y,y] are—though seldom encountered—quite tedious to type. Mathematica recognized this alternative:

Derivative[4,6][h][x,y]=== D[h[x,y],x,x,x,x,y,y,y,y,y,y]

Unspecified Functions

Suppose we write f[x] and g[x] to indicate that we are thinking about a couple of generic functions—functions to which we have assigned no specific meaning. Mathematica is content to work at such a level of abstraction, and is easily able to return all the standard relationships—identities—in, if you request, generalized form:

Clear[f,g]

D[f[x]g[x],x] (*Product Rule*)

D[f[x]g[x],{x,2}] (*Leibniz' Rule: generalized Product Rule*)
D[f[x]g[x],{x,3}]
D[f[x]g[x],{x,4}]

D[f[x]/g[x],x] (*Quotient Rule*)

Together[%]

D[f[g[x]],x] (*Chain Rule and higher order generalizations*)
D[f[g[x]],{x,2}] 
D[f[g[x]],{x,3}]

NOTE the notation that Mathematica reserves for derivatives of order higher than 2.

D[(∫_0^x)f[u]\[DifferentialD]u,x] (*Fundamental Theorem of the Calculus*)
(∫_0^x)D[f[u],u]\[DifferentialD]u

Mathematica is comfortable also with dummy functions of several variables.  For example, we have

D[f[x,y]g[x,y],x,y]
D[f[x,g[x,y]],x,y]

NOTE the notation used by Mathematica to name partial derivatives.

Implicit Differentiation

This technique is often quite useful, though it is my impression that it is no longer taught in introductory calculus courses.  I learned it from F. L. Griffin (Reed College's very first faculty member, in whose last class (Analytical Geometry, Spring 1952) I was a student), and take my example from §30 of his Mathematical Analysis Higher Course (1927).

Clear[x,y,a]

Suppose we know of y[x] only that (x^2)+(y^2)=10x y+(a^2) and want to obtain y''[x]. Rather than solving for y[x]and then double differentiating, we proceed by clever indirection: we construct

D[(x^2)+((y[x])^2)-10x y[x]-(a^2),x]

surely vanishes; we use this information to obtain a description of y'[x]:

Solve[2 x-10 y[x]-10 x y'[x]+2 y[x] y'[x]==0,y'[x]]

Define

firstderivative[x_]:=((x-5 y[x])/(5 x-y[x]))

and obtain

D[firstderivative[x],x]

Define

secondderivative[x_,p_]:=(((-x+5 y[x]) (5-p))/(((5 x-y[x]))^2))-((-1+5 p)/(5 x-y[x]))

and obtain

secondderivative[x,firstderivative[x]]

so0

Simplify[%]

But it was announced at the outset that (x^2)+(y^2)-10x y=(a^2) so we have finally

y''[x]==((24(a^2))/(((y-5x))^3))

Maxima & Minima

The following commands serve to define a simple function and to locate its extrema:

Clear[f]

f[x_]:=(x^3)-9x+5

Plot[f[x],{x,-4,4}]

Solve[f'[x]==0,x]

N[%]

Those results agree with the impression conveyed by the figure.  We will decorate the figure to make the point more vivid:

extrema={x,f[x]}/.%

Plot[f[x],{x,-4,4},
Epilog->{PointSize[0.03],Hue[1],Map[Point, extrema]}]

I have borrowed this discussion from Torrence & Torrence (The Students' Introduction to Mathematica (1999)), page 140.  Concerning the blue details, they have this to say: The Epilog command can be used with any command that produces graphics…It allows you to overlay "graphics primitives" (such as points) on the graphic after it has been rendered. [We adjusted the pointsize and color to make them more visible.] The Map[Point, extrema] produces a list of primitive point objects, one for each point in the list we named extrema.

Concerning the ReplaceAll command /. see

?/.

Limits & Series

These are distinct but closely related topics, which I have lumped  together for reasons that in Mathematica 6 have become irrelevant.  For basic information, command

?Limit
?Series

and consult the sites to which » leads you.

Evaluation of Limits

First Example

Clear[f]

f[x_]:=((Sin[x])/x)

At x = 0 we get the improper construction 0/0, yet from the figure

Plot[f[x],{x,-6π,6π}, PlotRange->All, Ticks->{{-6π,-4π,-2π,2π,4π,6π},{1}}]

we see that the function is perfectly well-behaved there, and are not surprised to obtain

Limit[f[x],x->0]

Expansion about the origin shows how this result is achieved:

Series[f[x], {x,0,4}]

Also of interest is the statement

Limit[f[x],x->∞]

Second Example

Clear[g]

Look now to the function

g[x_]:=(1/x)

Plot[g[x],{x,-1,1}, PlotRange->{-100,100}]

Limit[g[x],x->0]

Evidently Mathematica has, by default, elected to approach the origin from the right. One can reverse the direction of approach by modifying the command:

Limit[g[x],x->0, Direction->1]

Third Example

While the function

(E^(-x))

has well defined asymptotic limits

Limit[(E^(-x)),x->∞]
Limit[(E^(-x)),x->-∞]

the function

Sin[x]

possesses no such limits, but forever buzzes within a well-defined interval. Mathematica tells you so:

Limit[Sin[x],x->∞]
Limit[Sin[x],x->-∞]

Construction and Summation of Power Series

These are obviously related—yet quite distinct—subjects. To discover what Mathematica has to say about them, go to Documentation Center and search for tutorial/SeriesLimitsAndResidues/Overview.

Construction of Series

Examples: Some Familiar, Some Less Familiar

Series[Log[1+x],{x,0,10}]
Series[ArcTan[x],{x,0,10}]
Series[(1/(1-x)),{x,0,10}]
Series[(1/(1-x))Log[(1/(1-x))],{x,0,10}]

In the latter connection we note that

Table[Sum[(1/k),{k,1,n}],{n,1,10}]

Just as interesting integrals are sometimes taken to define numbers/functions, so are interesting series. For example, the following series serves to define the Bernoulli numbers

Series[(t/((E^t)-1)),{t,0,10}]

Table[((BernoulliB[n])/(n!)),{n,0,10}]

The series discussed above are all Maclaurin series (expansions about x = 0). Taylor series are produced by a slight variation of the same command:

Clear[f]

Series[f[x],{x,a,4}]

REMARK: The command

Normal[%]

lops off the annoying  ((O[ ])^( n)) term that otherwise terminates every application of the Series command.

Sums, and the Summation of Series

Examples: Some Familiar, Some Less Familiar

(∑_(k=1)^n)k
(∑_(k=1)^n)(k^2)
(∑_(k=1)^n)(k^3)

Just as definite integrals very often serve to define the so-called "higher functions," so also do infinite sums:

(∑_(n=0)^∞)(1/(n!))(x^n)
(∑_(n=0)^∞)(1/(((n!))^2))(x^n)
(∑_(n=0)^∞)(1/((2n)!))(x^n)
(∑_(n=0)^∞)((n!)/((2n)!))(x^n)
(∑_(n=1)^∞)(1/(n^x))

REMARK: The keyboard command EscsumtEsc serves to create

(∑_(\[Placeholder]=\[Placeholder])^\[Placeholder])\[Placeholder]